Fit VBGE models to back-calculated length-at-age

Author

Max Lindmark, Jan Ohlberger, Anna Gårdmark

Published

July 4, 2023

Load libraries

pkgs <- c("tidyverse", "tidylog", "rTPC", "nls.multstart", "broom", "RColorBrewer", 
          "viridis", "minpack.lm", "patchwork", "ggtext", "brms", "tidybayes", "modelr")

# minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library,character.only = T))
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'tidylog'


The following objects are masked from 'package:dplyr':

    add_count, add_tally, anti_join, count, distinct, distinct_all,
    distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
    full_join, group_by, group_by_all, group_by_at, group_by_if,
    inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
    relocate, rename, rename_all, rename_at, rename_if, rename_with,
    right_join, sample_frac, sample_n, select, select_all, select_at,
    select_if, semi_join, slice, slice_head, slice_max, slice_min,
    slice_sample, slice_tail, summarise, summarise_all, summarise_at,
    summarise_if, summarize, summarize_all, summarize_at, summarize_if,
    tally, top_frac, top_n, transmute, transmute_all, transmute_at,
    transmute_if, ungroup


The following objects are masked from 'package:tidyr':

    drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
    spread, uncount


The following object is masked from 'package:stats':

    filter


Loading required package: viridisLite

Loading required package: Rcpp

Loading 'brms' package (version 2.19.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').


Attaching package: 'brms'


The following object is masked from 'package:stats':

    ar



Attaching package: 'tidybayes'


The following objects are masked from 'package:brms':

    dstudent_t, pstudent_t, qstudent_t, rstudent_t



Attaching package: 'modelr'


The following object is masked from 'package:broom':

    bootstrap
# devtools::install_github("seananderson/ggsidekick") ## not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Load functions
home <- here::here()
fxn <- list.files(paste0(home, "/R/functions"))
invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))

Load cache

# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/analyze-data/02-fit-vbge_cache/html"))

Read and trim data

d <- #read_csv(paste0(home, "/data/for-analysis/dat.csv")) |> 
  read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv") |> 
  filter(age_ring == "Y", # use only length-at-age by filtering on age_ring
         !area == "VN") 
Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 36,279 rows (11%), 302,181 rows remaining
# Sample size by area and cohort
ns <- d |> 
  group_by(cohort, area) |> 
  summarise(n = n())
group_by: 2 grouping variables (cohort, area)
summarise: now 453 rows and 3 columns, one group variable remaining (cohort)
# Minimum number of observations per area and cohort
d$area_cohort <- as.factor(paste(d$area, d$cohort))

d <- d |>
  group_by(area_cohort) |> 
  filter(n() > 100)
group_by: one grouping variable (area_cohort)
filter (grouped): removed 2,796 rows (1%), 299,385 rows remaining
# Minimum number of observations per area, cohort, and age
d$area_cohort_age <- as.factor(paste(d$area, d$cohort, d$age_bc))

d <- d |>
  group_by(area_cohort_age) |> 
  filter(n() > 10)
group_by: one grouping variable (area_cohort_age)
filter (grouped): removed 3,584 rows (1%), 295,801 rows remaining
# Minimum number of cohorts in a given area
cnt <- d |>
  group_by(area) |>
  summarise(n = n_distinct(cohort)) |>
  filter(n >= 10)
group_by: one grouping variable (area)
summarise: now 11 rows and 2 columns, ungrouped
filter: no rows removed
d <- d[d$area %in% cnt$area, ]

# Plot cleaned data
ggplot(d, aes(age_bc, length_mm)) +
  geom_jitter(size = 0.1, height = 0, alpha = 0.1) +
  scale_x_continuous(breaks = seq(20)) +
  theme(axis.text.x = element_text(angle = 0)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  labs(x = "Age", y = "Length (mm)") +
  guides(color = "none") + 
  facet_wrap(~area, scale = "free_x")

# Longitude and latitude attributes for each area
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) |>
  mutate_at(c("lat", "lon"), as.numeric)
mutate_at: converted 'lat' from character to double (0 new NA)
           converted 'lon' from character to double (0 new NA)

Fit VBGE models

# Get individual growth parameters (functions: VBGF/Gompertz,nls_out,fit_nls)
IVBG <- d |> 
  group_by(ID) |> 
  summarize(k = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$k,
            k_se = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$k_se,
            linf = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$linf,
            linf_se = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$linf_se)
group_by: one grouping variable (ID)
summarize: now 79,122 rows and 5 columns, ungrouped

Inspect predictions

IVBG <- IVBG |> drop_na(k) # The NAs are because the model didn't converge or because they were below the threshold age
drop_na: removed 41,442 rows (52%), 37,680 rows remaining
IVBG <- IVBG |>
  mutate(k_lwr = k - 1.96*k_se,
         k_upr = k + 1.96*k_se,
         linf_lwr = linf - 1.96*linf_se,
         linf_upr = linf + 1.96*linf_se,
         row_id = row_number())
mutate: new variable 'k_lwr' (double) with 37,641 unique values and 0% NA
        new variable 'k_upr' (double) with 37,641 unique values and 0% NA
        new variable 'linf_lwr' (double) with 37,641 unique values and 0% NA
        new variable 'linf_upr' (double) with 37,641 unique values and 0% NA
        new variable 'row_id' (integer) with 37,680 unique values and 0% NA
# Plot all K's
IVBG |>
  #filter(row_id < 5000) |>
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr)) +
  geom_point(alpha = 0.5) +
  geom_errorbar(alpha = 0.5) +
  NULL

# Plot all L_inf's
IVBG |>
  #filter(row_id < 5000) |>
  ggplot(aes(row_id, linf, ymin = linf_lwr, ymax = linf_upr)) +
  geom_point(alpha = 0.5) +
  geom_errorbar(alpha = 0.5) +
  NULL

# We can also consider removing individuals where the SE of k is larger than the fit
IVBG |>
  #mutate(keep = ifelse(k > quantile(k_se, probs = 0.95), "N", "Y")) |>
  mutate(keep = ifelse(k < k_se, "N", "Y")) |>
  #filter(row_id < 10000) |>
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~keep, ncol = 1) +
  geom_errorbar(alpha = 0.5) +
  NULL
mutate: new variable 'keep' (character) with 2 unique values and 0% NA

# Add back cohort and area variables
IVBG <- IVBG |> 
  left_join(d |> select(ID, area_cohort)) |> 
  separate(area_cohort, into = c("area", "cohort"), remove = TRUE, sep = " ") |> 
  mutate(cohort = as.numeric(cohort))
Adding missing grouping variables: `area_cohort_age`
select: dropped 11 variables (length_mm, age_bc, age_catch, age_ring, area, …)
Joining with `by = join_by(ID)`
left_join: added 2 columns (area_cohort_age, area_cohort)
> rows only in x 0
> rows only in y ( 97,478)
> matched rows 198,323 (includes duplicates)
> =========
> rows total 198,323
mutate: converted 'cohort' from character to double (0 new NA)
# Summarise and save for sample size
sample_size <- IVBG |>
  group_by(area) |> 
  summarise(n_cohort = length(unique(cohort)),
            min_cohort = min(cohort),
            max_cohort = min(cohort),
            n_individuals = length(unique(ID)),
            n_data_points = n())
group_by: one grouping variable (area)
summarise: now 11 rows and 6 columns, ungrouped
sample_size
# A tibble: 11 × 6
   area  n_cohort min_cohort max_cohort n_individuals n_data_points
   <chr>    <int>      <dbl>      <dbl>         <int>         <int>
 1 BS          18       1985       1985          2370         11832
 2 BT          27       1972       1972          1882          9055
 3 FB          48       1969       1969          8678         47933
 4 FM          38       1963       1963          6704         39230
 5 HO          34       1982       1982          2711         12454
 6 JM          63       1953       1953          7738         40229
 7 MU          19       1981       1981          1920          9930
 8 RA          15       1985       1985          1205          5654
 9 SI_EK       35       1968       1968          1396          6523
10 SI_HA       41       1967       1967          2991         15143
11 TH           5       2000       2000            85           340
sample_size |>
  ungroup() |>
  summarise(sum_ind = sum(n_individuals), sum_n = sum(n_data_points))
ungroup: no grouping variables
summarise: now one row and 2 columns, ungrouped
# A tibble: 1 × 2
  sum_ind  sum_n
    <int>  <int>
1   37680 198323
write_csv(sample_size, paste0(home, "/output/sample_size.csv"))

# Compare how the means and quantiles differ depending on this filtering
IVBG_filter <- IVBG |>
  drop_na(k_se) |>
  #filter(k_se < quantile(k_se, probs = 0.95)) |> 
  filter(k_se < k)
drop_na: no rows removed
filter: removed 10,208 rows (5%), 188,115 rows remaining
# Summarize growth coefficients by cohort and area
VBG <- IVBG |>
  filter(k_se < k) |> # new!
  group_by(cohort, area) |>
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            linf_median = quantile(linf, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) |> 
  ungroup()
filter: removed 10,208 rows (5%), 188,115 rows remaining
group_by: 2 grouping variables (cohort, area)
summarize: now 343 rows and 7 columns, one group variable remaining (cohort)
ungroup: no grouping variables
VBG_filter <- IVBG_filter |>
  group_by(cohort, area) |>
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) |> 
  ungroup()
group_by: 2 grouping variables (cohort, area)
summarize: now 343 rows and 6 columns, one group variable remaining (cohort)
ungroup: no grouping variables
ggplot() +
  geom_ribbon(data = VBG, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,
                              fill = "All k's"), alpha = 0.4) +
  geom_ribbon(data = VBG_filter, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,
                                     fill = "Filtered"), alpha = 0.4) +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered")) + 
  guides(fill = "none") +
  facet_wrap(~area)

ggplot() +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered")) + 
  guides(fill = "none") +
  facet_wrap(~area)

# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.

BT 1996 has a very high k, so I’ll inspect it in more detail

bt96 <- IVBG |> 
  filter(area == "BT" & cohort %in% c(1995, 1996, 1997, 1998))
filter: removed 197,958 rows (>99%), 365 rows remaining
ggplot(bt96, aes(as.factor(cohort), k)) + 
  geom_jitter(width = 0.2, height = 0) + 
  geom_boxplot(fill = NA, width = 0.2) + 
  geom_hline(yintercept = mean(filter(bt96, cohort == 1996)$k), linetype = 2)
filter: removed 293 rows (80%), 72 rows remaining

# Check the fits in more detail
bt96b <- bt96 |> filter(cohort == 1996) |> as.data.frame()
filter: removed 293 rows (80%), 72 rows remaining
t <- bt96b |> 
  filter(k_se < k) |>
  filter(linf_se < linf) |> 
  filter(linf_lwr > 0)
filter: removed 16 rows (22%), 56 rows remaining
filter: no rows removed
filter: no rows removed
t
                        ID         k       k_se     linf    linf_se       k_lwr
1  2000_139_BT.female.4.NA 0.6998057 0.03282194 176.1500   3.041886  0.63547474
2  2000_139_BT.female.4.NA 0.6998057 0.03282194 176.1500   3.041886  0.63547474
3  2000_139_BT.female.4.NA 0.6998057 0.03282194 176.1500   3.041886  0.63547474
4  2000_139_BT.female.4.NA 0.6998057 0.03282194 176.1500   3.041886  0.63547474
5  2000_169_BT.female.4.NA 0.8309359 0.28906869 178.6083  20.099801  0.26436123
6  2000_169_BT.female.4.NA 0.8309359 0.28906869 178.6083  20.099801  0.26436123
7  2000_169_BT.female.4.NA 0.8309359 0.28906869 178.6083  20.099801  0.26436123
8  2000_169_BT.female.4.NA 0.8309359 0.28906869 178.6083  20.099801  0.26436123
9   2000_24_BT.female.4.NA 0.5728722 0.05217399 176.2271   6.819268  0.47061116
10  2000_24_BT.female.4.NA 0.5728722 0.05217399 176.2271   6.819268  0.47061116
11  2000_24_BT.female.4.NA 0.5728722 0.05217399 176.2271   6.819268  0.47061116
12  2000_24_BT.female.4.NA 0.5728722 0.05217399 176.2271   6.819268  0.47061116
13   2000_4_BT.female.4.NA 0.2457892 0.06377208 435.4096  75.698097  0.12079595
14   2000_4_BT.female.4.NA 0.2457892 0.06377208 435.4096  75.698097  0.12079595
15   2000_4_BT.female.4.NA 0.2457892 0.06377208 435.4096  75.698097  0.12079595
16   2000_4_BT.female.4.NA 0.2457892 0.06377208 435.4096  75.698097  0.12079595
17  2000_50_BT.female.4.NA 0.5000396 0.13927044 318.8768  41.323942  0.22706959
18  2000_50_BT.female.4.NA 0.5000396 0.13927044 318.8768  41.323942  0.22706959
19  2000_50_BT.female.4.NA 0.5000396 0.13927044 318.8768  41.323942  0.22706959
20  2000_50_BT.female.4.NA 0.5000396 0.13927044 318.8768  41.323942  0.22706959
21  2000_51_BT.female.4.NA 0.6280041 0.19388224 206.1642  25.353579  0.24799490
22  2000_51_BT.female.4.NA 0.6280041 0.19388224 206.1642  25.353579  0.24799490
23  2000_51_BT.female.4.NA 0.6280041 0.19388224 206.1642  25.353579  0.24799490
24  2000_51_BT.female.4.NA 0.6280041 0.19388224 206.1642  25.353579  0.24799490
25  2000_57_BT.female.4.NA 0.3710587 0.19779509 333.8080  98.764561 -0.01661973
26  2000_57_BT.female.4.NA 0.3710587 0.19779509 333.8080  98.764561 -0.01661973
27  2000_57_BT.female.4.NA 0.3710587 0.19779509 333.8080  98.764561 -0.01661973
28  2000_57_BT.female.4.NA 0.3710587 0.19779509 333.8080  98.764561 -0.01661973
29  2000_78_BT.female.4.NA 0.4577937 0.05308862 208.3834  11.887232  0.35374004
30  2000_78_BT.female.4.NA 0.4577937 0.05308862 208.3834  11.887232  0.35374004
31  2000_78_BT.female.4.NA 0.4577937 0.05308862 208.3834  11.887232  0.35374004
32  2000_78_BT.female.4.NA 0.4577937 0.05308862 208.3834  11.887232  0.35374004
33  2000_82_BT.female.4.NA 0.5529989 0.06526315 179.5668   9.224169  0.42508312
34  2000_82_BT.female.4.NA 0.5529989 0.06526315 179.5668   9.224169  0.42508312
35  2000_82_BT.female.4.NA 0.5529989 0.06526315 179.5668   9.224169  0.42508312
36  2000_82_BT.female.4.NA 0.5529989 0.06526315 179.5668   9.224169  0.42508312
37  2000_83_BT.female.4.NA 0.6063114 0.05190447 178.8868   6.254236  0.50457860
38  2000_83_BT.female.4.NA 0.6063114 0.05190447 178.8868   6.254236  0.50457860
39  2000_83_BT.female.4.NA 0.6063114 0.05190447 178.8868   6.254236  0.50457860
40  2000_83_BT.female.4.NA 0.6063114 0.05190447 178.8868   6.254236  0.50457860
41  2000_95_BT.female.4.NA 0.1261268 0.07828577 774.2764 390.036643 -0.02731335
42  2000_95_BT.female.4.NA 0.1261268 0.07828577 774.2764 390.036643 -0.02731335
43  2000_95_BT.female.4.NA 0.1261268 0.07828577 774.2764 390.036643 -0.02731335
44  2000_95_BT.female.4.NA 0.1261268 0.07828577 774.2764 390.036643 -0.02731335
45 2001_201_BT.female.5.NA 0.5572245 0.11202672 366.7864  31.931152  0.33765214
46 2001_201_BT.female.5.NA 0.5572245 0.11202672 366.7864  31.931152  0.33765214
47 2001_201_BT.female.5.NA 0.5572245 0.11202672 366.7864  31.931152  0.33765214
48 2001_201_BT.female.5.NA 0.5572245 0.11202672 366.7864  31.931152  0.33765214
49 2001_222_BT.female.5.NA 0.4204646 0.14090368 397.1261  68.884483  0.14429342
50 2001_222_BT.female.5.NA 0.4204646 0.14090368 397.1261  68.884483  0.14429342
51 2001_222_BT.female.5.NA 0.4204646 0.14090368 397.1261  68.884483  0.14429342
52 2001_222_BT.female.5.NA 0.4204646 0.14090368 397.1261  68.884483  0.14429342
53 2001_236_BT.female.5.NA 0.2093150 0.09910580 573.1898 192.597276  0.01506767
54 2001_236_BT.female.5.NA 0.2093150 0.09910580 573.1898 192.597276  0.01506767
55 2001_236_BT.female.5.NA 0.2093150 0.09910580 573.1898 192.597276  0.01506767
56 2001_236_BT.female.5.NA 0.2093150 0.09910580 573.1898 192.597276  0.01506767
       k_upr   linf_lwr  linf_upr row_id area_cohort_age area cohort
1  0.7641367 170.187889  182.1121  25583       BT 1996 1   BT   1996
2  0.7641367 170.187889  182.1121  25583       BT 1996 2   BT   1996
3  0.7641367 170.187889  182.1121  25583       BT 1996 3   BT   1996
4  0.7641367 170.187889  182.1121  25583       BT 1996 4   BT   1996
5  1.3975105 139.212721  218.0039  25667       BT 1996 1   BT   1996
6  1.3975105 139.212721  218.0039  25667       BT 1996 2   BT   1996
7  1.3975105 139.212721  218.0039  25667       BT 1996 3   BT   1996
8  1.3975105 139.212721  218.0039  25667       BT 1996 4   BT   1996
9  0.6751332 162.861337  189.5929  25905       BT 1996 1   BT   1996
10 0.6751332 162.861337  189.5929  25905       BT 1996 2   BT   1996
11 0.6751332 162.861337  189.5929  25905       BT 1996 3   BT   1996
12 0.6751332 162.861337  189.5929  25905       BT 1996 4   BT   1996
13 0.3707825 287.041290  583.7778  26101       BT 1996 1   BT   1996
14 0.3707825 287.041290  583.7778  26101       BT 1996 2   BT   1996
15 0.3707825 287.041290  583.7778  26101       BT 1996 3   BT   1996
16 0.3707825 287.041290  583.7778  26101       BT 1996 4   BT   1996
17 0.7730097 237.881859  399.8717  26103       BT 1996 1   BT   1996
18 0.7730097 237.881859  399.8717  26103       BT 1996 2   BT   1996
19 0.7730097 237.881859  399.8717  26103       BT 1996 3   BT   1996
20 0.7730097 237.881859  399.8717  26103       BT 1996 4   BT   1996
21 1.0080133 156.471211  255.8572  26104       BT 1996 1   BT   1996
22 1.0080133 156.471211  255.8572  26104       BT 1996 2   BT   1996
23 1.0080133 156.471211  255.8572  26104       BT 1996 3   BT   1996
24 1.0080133 156.471211  255.8572  26104       BT 1996 4   BT   1996
25 0.7587370 140.229501  527.3866  26119       BT 1996 1   BT   1996
26 0.7587370 140.229501  527.3866  26119       BT 1996 2   BT   1996
27 0.7587370 140.229501  527.3866  26119       BT 1996 3   BT   1996
28 0.7587370 140.229501  527.3866  26119       BT 1996 4   BT   1996
29 0.5618474 185.084375  231.6823  26157       BT 1996 1   BT   1996
30 0.5618474 185.084375  231.6823  26157       BT 1996 2   BT   1996
31 0.5618474 185.084375  231.6823  26157       BT 1996 3   BT   1996
32 0.5618474 185.084375  231.6823  26157       BT 1996 4   BT   1996
33 0.6809147 161.487402  197.6461  26163       BT 1996 1   BT   1996
34 0.6809147 161.487402  197.6461  26163       BT 1996 2   BT   1996
35 0.6809147 161.487402  197.6461  26163       BT 1996 3   BT   1996
36 0.6809147 161.487402  197.6461  26163       BT 1996 4   BT   1996
37 0.7080441 166.628541  191.1451  26166       BT 1996 1   BT   1996
38 0.7080441 166.628541  191.1451  26166       BT 1996 2   BT   1996
39 0.7080441 166.628541  191.1451  26166       BT 1996 3   BT   1996
40 0.7080441 166.628541  191.1451  26166       BT 1996 4   BT   1996
41 0.2795669   9.804558 1538.7482  26194       BT 1996 1   BT   1996
42 0.2795669   9.804558 1538.7482  26194       BT 1996 2   BT   1996
43 0.2795669   9.804558 1538.7482  26194       BT 1996 3   BT   1996
44 0.2795669   9.804558 1538.7482  26194       BT 1996 4   BT   1996
45 0.7767969 304.201370  429.3715  26927       BT 1996 1   BT   1996
46 0.7767969 304.201370  429.3715  26927       BT 1996 2   BT   1996
47 0.7767969 304.201370  429.3715  26927       BT 1996 3   BT   1996
48 0.7767969 304.201370  429.3715  26927       BT 1996 4   BT   1996
49 0.6966359 262.112502  532.1397  27021       BT 1996 1   BT   1996
50 0.6966359 262.112502  532.1397  27021       BT 1996 2   BT   1996
51 0.6966359 262.112502  532.1397  27021       BT 1996 3   BT   1996
52 0.6966359 262.112502  532.1397  27021       BT 1996 4   BT   1996
53 0.4035624 195.699170  950.6805  27079       BT 1996 1   BT   1996
54 0.4035624 195.699170  950.6805  27079       BT 1996 2   BT   1996
55 0.4035624 195.699170  950.6805  27079       BT 1996 3   BT   1996
56 0.4035624 195.699170  950.6805  27079       BT 1996 4   BT   1996
bt96 |> 
  filter(k_se < k) |>
  filter(linf_se < linf) |> 
  filter(linf_lwr > 0) |> 
  ggplot(aes(as.factor(cohort), k)) + 
  geom_jitter(width = 0.2, height = 0) + 
  geom_boxplot(fill = NA, width = 0.2) + 
  geom_hline(yintercept = mean(filter(bt96, cohort == 1996)$k), linetype = 2)
filter: removed 45 rows (12%), 320 rows remaining
filter: no rows removed
filter: removed 33 rows (10%), 287 rows remaining
filter: removed 293 rows (80%), 72 rows remaining

# Plot these individuals in t
t <- t |> 
  left_join(d, by = "ID") |> 
  mutate(length_mm_pred = linf*(1-exp(-k*age_bc)))
left_join: added 16 columns (area_cohort_age.x, area.x, cohort.x, length_mm, age_bc, …)
           > rows only in x         0
           > rows only in y  (295,745)
           > matched rows         224    (includes duplicates)
           >                 =========
           > rows total           224
mutate: new variable 'length_mm_pred' (double) with 56 unique values and 0% NA
ggplot(t, aes(age_bc, length_mm)) +
  geom_jitter(width = 0.3) +
  geom_line(aes(age_bc, length_mm_pred)) +
  facet_wrap(~ID)

# These fits look totally reasonable

Add GAM-predicted temperature to growth models

pred_temp <- read_csv(paste0(home, "/output/gam_predicted_temps.csv")) |> 
  rename(cohort = year)
Rows: 395 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): area, model
dbl (2): year, temp

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (cohort)
VBG <- VBG |>
  left_join(pred_temp, by = c("area", "cohort"))
left_join: added 2 columns (temp, model)
           > rows only in x     0
           > rows only in y  ( 52)
           > matched rows     343
           >                 =====
           > rows total       343
# Save data for map-plot
cohort_sample_size <- IVBG |>
  group_by(area, cohort) |> 
  summarise(n = n()) # individuals, not samples!
group_by: 2 grouping variables (area, cohort)
summarise: now 343 rows and 3 columns, one group variable remaining (area)
VBG <- left_join(VBG, cohort_sample_size, by = c("cohort", "area"))
left_join: added one column (n)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     343
           >                 =====
           > rows total       343
write_csv(VBG, paste0(home, "/output/vbg.csv"))

# Calculate the plotting order (also used for map plot)
order <- VBG |>
  ungroup() |>
  group_by(area) |>
  summarise(min_temp = min(temp)) |>
  arrange(desc(min_temp))
ungroup: no grouping variables
group_by: one grouping variable (area)
summarise: now 11 rows and 2 columns, ungrouped
order
# A tibble: 11 × 2
   area  min_temp
   <chr>    <dbl>
 1 SI_HA     7.83
 2 TH        6.85
 3 BS        6.22
 4 BT        5.87
 5 FM        5.87
 6 SI_EK     5.49
 7 MU        5.16
 8 FB        5.04
 9 HO        3.99
10 JM        3.37
11 RA        3.06
nareas <- length(unique(order$area)) + 2 # to skip the brightest colors
colors <- colorRampPalette(brewer.pal(name = "RdYlBu", n = 10))(nareas)[-c(7,8)]

Plot VBGE fits

# Sample 50 IDs per area and plot their data and VBGE fits
set.seed(4)
ids <- IVBG |> distinct(ID, .keep_all = TRUE) |> group_by(area) |> slice_sample(n = 30)
distinct: removed 160,643 rows (81%), 37,680 rows remaining
group_by: one grouping variable (area)
slice_sample (grouped): removed 37,350 rows (99%), 330 rows remaining
#ids |> ungroup() |> group_by(area) |> summarise(n = length(unique(ID))) |> arrange(n)

IVBG2 <- IVBG |>
  filter(ID %in% ids$ID) |> 
  distinct(ID, .keep_all = TRUE) |> 
  select(ID, k, linf)
filter: removed 196,691 rows (99%), 1,632 rows remaining
distinct: removed 1,302 rows (80%), 330 rows remaining
select: dropped 10 variables (k_se, linf_se, k_lwr, k_upr, linf_lwr, …)
d2 <- d |>
  ungroup() |>
  filter(ID %in% ids$ID) |>
  left_join(IVBG2, by = "ID") |>
  mutate(length_mm_pred = linf*(1-exp(-k*age_bc)))
ungroup: no grouping variables
filter: removed 294,169 rows (99%), 1,632 rows remaining
left_join: added 2 columns (k, linf)
           > rows only in x       0
           > rows only in y  (    0)
           > matched rows     1,632
           >                 =======
           > rows total       1,632
mutate: new variable 'length_mm_pred' (double) with 1,632 unique values and 0% NA
ggplot(d2, aes(age_bc, length_mm, group = ID, color = ID)) +
  geom_jitter(width = 0.3, height = 0, alpha = 0.5, size = 0.4) +
  geom_line(aes(age_bc, length_mm_pred, group = ID, color = ID),
            inherit.aes = FALSE, alpha = 0.8, linewidth = 0.3) + 
  guides(color = "none") +
  scale_color_viridis(discrete = TRUE, name = "Area", option = "cividis") +
  labs(x = "Age (yrs)", y = "Length (mm)") +
  facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area), ncol = 3)

ggsave(paste0(home, "/figures/vb_fits.pdf" ), width = 17, height = 17, units = "cm")

# rugs <- IVBG |> 
#   group_by(area) |> 
#   summarise(median_k = median(k),
#             median_linf = median(linf))

# k <- IVBG |> 
#   ggplot(aes(k, color = factor(area, order$area))) + 
#   geom_density(alpha = 0.4, fill = NA, adjust = 1.5) +  
#   scale_color_manual(values = colors, name = "Area") +
#   coord_cartesian(expand = 0) + 
#   labs(x = expression(italic(k))) +
#   geom_rug(data = rugs, aes(median_k)) +
#   theme(legend.position = c(0.9, 0.5))

# violons instead?
k2 <- IVBG |> 
  ggplot(aes(factor(area, order$area), k, fill = factor(area, order$area))) + 
  geom_violin(alpha = 0.6, color = NA) +  
  geom_boxplot(outlier.color = NA, width = 0.2, alpha = 0.8, fill = NA, size = 0.4) +
  # geom_jitter(data = IVBG |> group_by(area) |> slice_sample(n = 100),
  #             aes(factor(area, order$area), k), height = 0, width = 0.2, alpha = 0.1,
  #             color = "grey40") +
  scale_fill_manual(values = colors, name = "Area") +
  guides(fill = "none", color = "none") +
  labs(x = "", y = expression(italic(k)))

k2

# linf <- IVBG |> 
#   ggplot(aes(linf, color = factor(area, order$area))) + 
#   geom_density(alpha = 0.4, fill = NA, adjust = 1.5) +  
#   scale_color_manual(values = colors, name = "Area") +
#   coord_cartesian(expand = 0, xlim = c(0, 2000)) + 
#   labs(x = expression(paste(italic(L[infinity]), " [cm]"))) +
#   geom_rug(data = rugs, aes(median_linf)) +
#   guides(color = "none")

linf2 <- IVBG |> 
  # group_by(area) |> 
  # mutate(upr = quantile(linf, probs = 0.95)) |> 
  # ungroup() |> 
  # filter(linf < upr) |> 
  ggplot(aes(factor(area, order$area), linf, fill = factor(area, order$area))) + 
  geom_violin(alpha = 0.6, color = NA) +  
  geom_boxplot(outlier.color = NA, width = 0.1, alpha = 0.8, fill = NA, size = 0.4) +
  coord_cartesian(ylim = c(0, 2300)) +
  # geom_jitter(data = IVBG |> group_by(area) |> slice_sample(n = 100),
  #             aes(factor(area, order$area), linf), height = 0, width = 0.2, alpha = 0.1,
  #             color = "grey40") +
  scale_fill_manual(values = colors, name = "Area") +
  guides(fill = "none", color = "none") +
  labs(x = "", y = expression(paste(italic(L[infinity]), " [cm]")))

linf2

#k / linf
k2 / linf2

ggsave(paste0(home, "/figures/vb_pars.pdf" ), width = 17, height = 17, units = "cm")

Fit Sharpe-Schoolfield model to K

By area

model <- 'sharpeschoolhigh_1981'

# Get starting values on full dataset for Sharpe-Schoolfield model
dat <- VBG |>
  select(k_median, temp, area) |>
  rename(rate = k_median)
select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
rename: renamed one variable (rate)
lower <- get_lower_lims(dat$temp, dat$rate, model_name = model)
upper <- get_upper_lims(dat$temp, dat$rate, model_name = model)
start <- get_start_vals(dat$temp, dat$rate, model_name = model)
  
# Sharpe-Schoolfield model fit to data for each area
preds <- NULL
pred <- NULL

for(a in unique(dat$area)) {
  
  # Get data
  dd <- dat |> filter(area == a)
  
  # Fit model
  fit <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dd,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )
  
  # Make predictions on new data
  new_data <- data.frame(temp = seq(min(dd$temp), max(dd$temp), length.out = 100))
  
  pred <- augment(fit, newdata = new_data) |> mutate(area = a)
  
  # Add to general data frame
  preds <- data.frame(rbind(preds, pred))
  
}
filter: removed 280 rows (82%), 63 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 305 rows (89%), 38 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 302 rows (88%), 41 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 308 rows (90%), 35 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 295 rows (86%), 48 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 316 rows (92%), 27 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 324 rows (94%), 19 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 309 rows (90%), 34 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 325 rows (95%), 18 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 328 rows (96%), 15 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA
filter: removed 338 rows (99%), 5 rows remaining
mutate: new variable 'area' (character) with one unique value and 0% NA

All areas pooled

# One sharpe schoolfield fitted to all data
fit_all <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dat,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )

# Make predictions on new data
new_data_all <- data.frame(temp = seq(min(dat$temp), max(dat$temp), length.out = 100))

pred_all <- augment(fit_all, newdata = new_data_all) |> 
  mutate(area = "all")
mutate: new variable 'area' (character) with one unique value and 0% NA
# Add t_opt
kb <- 8.62e-05
data.frame(par = names(coef(fit_all)), est = coef(fit_all)) |> 
  pivot_wider(names_from = par, values_from = est) |> 
  summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))
pivot_wider: reorganized (par, est) into (r_tref, e, eh, th) [was 4x2, now 1x4]
summarise: now one row and one column, ungrouped
# A tibble: 1 × 1
  t_opt
  <dbl>
1  6.53

Plot Sharpe Schoolfield fits

# Plot growth coefficients by year and area against mean SST

p1 <- preds |>
  ggplot(aes(temp, .fitted, color = factor(area, order$area))) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 1, alpha = 0.8) +
  geom_line(linewidth = 1) +
  geom_line(data = pred_all, aes(temp, .fitted), linewidth = 1, inherit.aes = FALSE, linetype = 2) +
  scale_color_manual(values = colors, name = "Area") +
  guides(color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  scale_x_continuous(breaks = seq(-5, 20, 1)) +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)") +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal")

p1 + facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area))

p1

ggsave(paste0(home, "/figures/sharpe_school.pdf" ), width = 17, height = 17, units = "cm")

Can we fit a single Sharpe Schoolfield with area-specific parameters with brms?

# Again, here are the data we are fitting:
ggplot(dat, aes(temp, rate, color = factor(area, order$area))) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6)  

# Which family to use? Well we can try with a Gaussian or a student-t first... 
# After that, maybe gamma or lognormal, if we predict negative K's
hist(dat$rate)

# Here's the equation:
sharpeschoolhigh_1981
function (temp, r_tref, e, eh, th, tref) 
{
    tref <- 273.15 + tref
    k <- 8.62e-05
    boltzmann.term <- r_tref * exp(e/k * (1/tref - 1/(temp + 
        273.15)))
    inactivation.term <- 1/(1 + exp(eh/k * (1/(th + 273.15) - 
        1/(temp + 273.15))))
    return(boltzmann.term * inactivation.term)
}
<bytecode: 0x7ff30d3cd380>
<environment: namespace:rTPC>
# Add in fixed parameters
dat$bk <- 8.62e-05
dat$tref <- 8 + 273.15

# We can use the nls() parameters as starting values
summary(fit_all)

Formula: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, 
    tref = 8)

Parameters:
       Estimate Std. Error t value Pr(>|t|)  
r_tref   0.5265     0.5818   0.905   0.3661  
e        1.3179     1.1940   1.104   0.2705  
eh       1.8364     0.7329   2.506   0.0127 *
th       6.5328     6.1852   1.056   0.2916  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05265 on 339 degrees of freedom

Number of iterations to convergence: 33 
Achieved convergence tolerance: 1.49e-08
# As a first guess, I use the nls estimate as the mean, and standard deviation something close to it. We can also bound som fo them by 0, since negative values shouldn't be possible
n=10000
hist(rnorm(0.4, 0.4, n=n), main = "rtref")

hist(rnorm(1.4, 0.6, n=n), main = "e") 

hist(rnorm(2.1, 1.5, n=n), main = "eh")

hist(rnorm(7.5, 4, n=n), main = "hh")

fit_prior <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1, 
     nl = TRUE),
  data = dat,
  cores = 2,
  chains = 2,
  iter = 1500,
  sample_prior = "only",
  seed = 9,
  prior = c(
    prior(normal(0.4, 0.4), nlpar = "rtref", lb = 0),
    prior(normal(1.4, 0.6), nlpar = "e"),
    prior(normal(2.1, 1.5), nlpar = "eh"),
    prior(normal(7.5, 4), nlpar = "th", lb = 0)
  )
)
Compiling Stan program...
Start sampling
# Global prior predictive check in relation to data. Doesn't loo too informative...
dat |>
  data_grid(temp = seq_range(temp, n = 51)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fit_prior) |> 
  mutate(rtref = 0.27, e = 0.7, eh = 1.8, th = 11.8) |> 
  mutate(man_pred = (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))) |> 
  ggplot(aes(temp, y = .epred)) +
  stat_lineribbon(.width = c(0.95), alpha = 0.3, color = "black", fill = "black") + 
  geom_line(aes(y = man_pred), color = "tomato3", linetype = 2) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area") +
  NULL
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
        new variable 'tref' (double) with one unique value and 0% NA
mutate (grouped): new variable 'rtref' (double) with one unique value and 0% NA
                  new variable 'e' (double) with one unique value and 0% NA
                  new variable 'eh' (double) with one unique value and 0% NA
                  new variable 'th' (double) with one unique value and 0% NA
mutate (grouped): new variable 'man_pred' (double) with 51 unique values and 0% NA

# Now make sure it converges with real data but no random effects
fit_global <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1,
     nl = TRUE),
  data = dat,
  cores = 4,
  chains = 4,
  iter = 4000,
  seed = 9,
  prior = c(
    prior(normal(0.4, 0.4), nlpar = "rtref", lb = 0),
    prior(normal(1.4, 0.6), nlpar = "e"),
    prior(normal(2.1, 1.5), nlpar = "eh"),
    prior(normal(7.5, 4), nlpar = "th", lb = 0)
  )
)
Compiling Stan program...
recompiling to avoid crashing R session
Start sampling
plot(fit_global)

pp_check(fit_global)
Using 10 posterior draws for ppc type 'dens_overlay' by default.

plot(conditional_effects(fit_global, effect = "temp"), points = TRUE)

dat |>
  data_grid(temp = seq_range(temp, n = 51)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fit_global) |> 
  mutate(rtref = 0.4391, e = 1.4258, eh = 2.1066, th = 7.5369) |> # from summary(fit_all)
  mutate(man_pred = (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))) |> 
  ggplot(aes(temp, y = .epred)) +
  stat_lineribbon(.width = c(0.95), alpha = 0.3, color = "black", fill = "black") + 
  geom_line(aes(y = man_pred), color = "tomato3", linetype = 2) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area")
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
        new variable 'tref' (double) with one unique value and 0% NA
mutate (grouped): new variable 'rtref' (double) with one unique value and 0% NA
                  new variable 'e' (double) with one unique value and 0% NA
                  new variable 'eh' (double) with one unique value and 0% NA
                  new variable 'th' (double) with one unique value and 0% NA
mutate (grouped): new variable 'man_pred' (double) with 51 unique values and 0% NA

Ok, seems to work with priors and everything. They are roughly as broad as can be and still have a fit. Now fit the full model, with random area effects, more iterations and chains!

# Still good enough, might change some priors. Now look at random area effects!
# Gaussian model
fitb <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1 + (1|area),
     nl = TRUE),
  data = dat,
  cores = 4,
  chains = 4,
  iter = 4000,
  sample_prior = "yes",
  save_pars = save_pars(all = TRUE),
  seed = 9,
  prior = c(
    prior(normal(0.4, 0.4), nlpar = "rtref", lb = 0),
    prior(normal(1.4, 0.6), nlpar = "e"),
    prior(normal(2.1, 1.5), nlpar = "eh"),
    prior(normal(7.5, 4), nlpar = "th", lb = 0)
  ),
  control = list(adapt_delta = 0.95, max_treedepth = 12)
)
Compiling Stan program...
Start sampling
# Student-t model
fitbs <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1 + (1|area),
     nl = TRUE),
  data = dat,
  family = student,
  cores = 4,
  chains = 4,
  iter = 4000,
  sample_prior = "yes",
  save_pars = save_pars(all = TRUE),
  seed = 9,
  prior = c(
    prior(normal(0.4, 0.4), nlpar = "rtref", lb = 0),
    prior(normal(1.4, 0.6), nlpar = "e"),
    prior(normal(2.1, 1.5), nlpar = "eh"),
    prior(normal(7.5, 4), nlpar = "th", lb = 0)
  ),
  control = list(adapt_delta = 0.95, max_treedepth = 12)
)
Compiling Stan program...
Start sampling

Compare fitted models based on ELPD. Preferred model in the first row!

fitb_loo <- loo(fitb, moment_match = TRUE)
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Warning: Found 2 observations with a pareto_k > 0.7 in model 'fitb'. It is
recommended to set 'reloo = TRUE' in order to calculate the ELPD without the
assumption that these observations are negligible. This will refit the model 2
times to compute the ELPDs for the problematic observations directly.
fitbs_loo <- loo(fitbs, moment_match = TRUE)
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
Warning: Found 1 observations with a pareto_k > 0.7 in model 'fitbs'. It is
recommended to set 'reloo = TRUE' in order to calculate the ELPD without the
assumption that these observations are negligible. This will refit the model 1
times to compute the ELPDs for the problematic observations directly.
loo_compare(fitb_loo, fitbs_loo)
      elpd_diff se_diff
fitbs   0.0       0.0  
fitb  -10.9      10.8

The student model is preferred. Inspect it!

# Check fit
summary(fitbs)
 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 343) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Group-Level Effects: 
~area (Number of levels: 11) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(rtref_Intercept)     0.06      0.04     0.00     0.16 1.00     1640     2267
sd(e_Intercept)         0.21      0.14     0.01     0.53 1.00     2193     2920
sd(eh_Intercept)        0.22      0.20     0.01     0.77 1.00     1768     2734
sd(th_Intercept)        0.72      0.62     0.03     2.30 1.00     1820     3322

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept     0.58      0.21     0.29     1.10 1.00     1705     2567
e_Intercept         1.14      0.41     0.42     1.97 1.00     1899     2678
eh_Intercept        1.63      0.44     0.71     2.42 1.00     2703     2948
th_Intercept        6.54      2.60     2.40    12.92 1.00     1666     2308

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.00     0.03     0.04 1.00     7536     6447
nu       10.49      4.47     5.19    21.66 1.00     8058     6057

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fitbs)

pp_check(fitbs)
Using 10 posterior draws for ppc type 'dens_overlay' by default.

ggsave(paste0(home, "/figures/supp/sharpe_brms_ppcheck.pdf" ), width = 11, height = 11, units = "cm")

# Plot prior vs posterior
post <- fitbs |>
  as_draws_df() |>
  dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) |> 
  rename(rtref = b_rtref_Intercept, 
         e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Posterior")
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (rtref, e, eh, th)
pivot_longer: reorganized (rtref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
prior <- fitbs |>
  as_draws_df() |>
  dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) |> 
    rename(rtref = prior_b_rtref, 
         e = prior_b_e,
         eh = prior_b_eh,
         th = prior_b_th) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Prior")
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 4 variables (rtref, e, eh, th)
pivot_longer: reorganized (rtref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]
mutate: new variable 'type' (character) with one unique value and 0% NA
dist <- bind_rows(prior, post)

ggplot(dist, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free") + 
  theme(legend.position = c(0.35, 0.9)) +
  labs(fill = "")

ggsave(paste0(home, "/figures/supp/sharpe_brms_prior.pdf" ), width = 17, height = 17, units = "cm")

# Calculate T_opt
fitbs |>
  as_draws_df() |>
  dplyr::select( b_e_Intercept, b_eh_Intercept, b_th_Intercept) |> 
  rename(e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) |> 
  mutate(kb = 8.62e-05) |> 
  ungroup() |> 
  mutate(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1))) |> 
  drop_na(t_opt) |> 
  summarise(lwr_t_opt = quantile(t_opt, probs = 0.1),
            median_t_opt = median(t_opt),
            upr_t_opt = quantile(t_opt, probs = 0.9))
Warning: Dropping 'draws_df' class as required metadata was removed.
rename: renamed 3 variables (e, eh, th)
mutate: new variable 'kb' (double) with one unique value and 0% NA
ungroup: no grouping variables
Warning: There was 1 warning in `.fun()`.
ℹ In argument: `t_opt = (eh * th)/(eh + kb * th * log((eh/e) - 1))`.
Caused by warning in `log()`:
! NaNs produced
mutate: new variable 't_opt' (double) with 7,871 unique values and 2% NA
drop_na: removed 130 rows (2%), 7,870 rows remaining
summarise: now one row and 3 columns, ungrouped
# A tibble: 1 × 3
  lwr_t_opt median_t_opt upr_t_opt
      <dbl>        <dbl>     <dbl>
1      3.76         6.10      10.1
# Make the main plot (conditional effect of temperature, with and without random effects)
# Predictions without random effects
glob <- dat |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs, re_formula = NA) |> 
  ungroup()
mutate: new variable 'bk' (double) with one unique value and 0% NA
        new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables
dat |>
  group_by(area) |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs) |> 
  ungroup() |> 
  ggplot(aes(temp, y = .epred, color = factor(area, order$area))) +
  stat_lineribbon(data = glob, aes(temp, .epred), .width = c(0.9), inherit.aes = FALSE,
                  fill = "black", color = "black", alpha = 0.1) +
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 1, alpha = 0.8) +
  stat_lineribbon(.width = c(0)) +
  stat_lineribbon(data = glob, aes(temp, .epred), .width = c(0), inherit.aes = FALSE,
                  color = "black", alpha = 0.9, linetype = 2) +
  scale_color_manual(values = colors, name = "Area") +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)")
group_by: one grouping variable (area)
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
        new variable 'tref' (double) with one unique value and 0% NA
ungroup: no grouping variables

ggsave(paste0(home, "/figures/sharpe_school_bayes.pdf" ), width = 17, height = 17, units = "cm")

Now plot comparisons with the nls multstart fits

# Compare with area-specific sharpe scool!
area_pred_brms <- dat |>
  group_by(area) |> 
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs)
group_by: one grouping variable (area)
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
        new variable 'tref' (double) with one unique value and 0% NA
# Plot area specific in comparison with nls multstart
p1 <- ggplot(area_pred_brms, aes(temp, y = .epred, color = factor(area, order$area),
                           fill = factor(area, order$area))) +
  stat_lineribbon(.width = c(0.95), alpha = 0.4, color = NA) +
  stat_lineribbon(.width = c(0)) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area") +
  scale_fill_manual(values = colors, name = "Area") +
  scale_linetype_manual(values = 2) +
  facet_wrap(~factor(area, rev(order$area))) +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       linetype = "")

p1

ggsave(paste0(home, "/figures/supp/sharpe_school_bayes_ci_facet.pdf" ), width = 17, height = 17, units = "cm")

# Now compare to nls fits
p1 + geom_line(data = preds, aes(temp, .fitted, linetype = "nls multstart"), color = "gray10")

# And for the population-level prediction
nd_brms <- dat |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15)
ungroup: no grouping variables
mutate: new variable 'bk' (double) with one unique value and 0% NA
        new variable 'tref' (double) with one unique value and 0% NA
pred_brms <- predict(fitbs, newdata = nd_brms, re_formula = NA) |> as.data.frame()
nd_brms$pred_pop <- pred_brms$Estimate

# (why not also comparing the global brms model...)
pred_brms_global <- predict(fit_global, newdata = nd_brms) |> as.data.frame()
nd_brms$pred_glob <- pred_brms_global$Estimate

ggplot(nd_brms, aes(temp, y = pred_pop, linetype = "mixed brms model")) +
  geom_line() + 
  geom_line(data = pred_all, aes(temp, .fitted, linetype = "nls"), linewidth = 1, inherit.aes = FALSE) +
  geom_line(data = nd_brms, aes(temp, pred_glob, linetype = "global brms model"), linewidth = 1, inherit.aes = FALSE) +
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area") +
  scale_linetype_manual(values = c(1, 2, 3)) + 
  guides(fill = "none",
         color = guide_legend(ncol = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1)),
         linetype = guide_legend(keywidth = unit(1, "cm"))) +
  theme(axis.title.y = ggtext::element_markdown()) +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       linetype = "Model")

Extra analysis

knitr::knit_exit()